!--------------------------------------------------!
!		MODULE MATRIXMOD		   !
!--------------------------------------------------!
! This module permits to operate with matrix       !
! in a direct way 				   !
!--------------------------------------------------!
! Written by Christian F. A. Negre, Csrdoba, UNC   !
! october 2007					   !
!--------------------------------------------------!  

module matrixmod

  use control
  use constants
  
  implicit none

  public :: Id,write_matrix,adj,product
  public :: diag, comm
  
  interface operator(.x.)
    module procedure product
  end interface

 contains
 
!--------------------------------------------------!
!	Identity Matrix				   !
!--------------------------------------------------!

  complex(dp) function Id()
    dimension :: Id(N,N)
    integer :: i
    Id=0.0_dp 
    do i=1,N
    Id(i,i)=1.0_dp
    enddo
  end function Id

!--------------------------------------------------!
!	Traspose conjugated Matrix		   !
!--------------------------------------------------!

  complex(dp) function adj(A)
    complex(dp) :: A(N,N)
    dimension :: adj(N,N)
    integer :: i,j
    do i=1,N
       do j=1,N
          adj(j,i) = conjg(A(i,j))
       enddo
    enddo
  end function adj

!--------------------------------------------------!
!	Writing a Matrix			   ! 
!--------------------------------------------------!

  subroutine write_matrix(A)
    complex(dp) :: A(N,N)
    integer :: i,j
    
    write(*,*)'***************************************************************'    
    write(*,*)'Real part'
    do i=1,N
      write(*,'(500f10.5)')(real(A(i,j)),j=1,N)	
    enddo
    write(*,*)'Imaginary part'
    do i=1,N
      write(*,'(500f10.5)')(aimag(A(i,j)),j=1,N)	
    enddo
    write(*,*)'***************************************************************'    
  end subroutine write_matrix  

!--------------------------------------------------!
!	Matrix Product				   !
!--------------------------------------------------!

  complex(dp) function product(A,B)
    complex(dp), intent(in) :: A(N,N), B(N,N)
    complex(dp) :: alpha, beta
!    integer  :: i,j,k
    dimension :: product(N,N)

!    do i=1,N
!       do j=1,N
!          product(i,j)= (0.0_dp,0.0_dp)
!          do k=1,N
!	     product(i,j) = product(i,j) + A(i,k)*B(k,j)
!          enddo
!       enddo
!    enddo

   alpha=1.0_dp
   beta=0.0_dp
   
   call zgemm('N', 'N', N, N, N, alpha, A, N, B, N, beta, product, N)  
     
  end function product

!--------------------------------------------------!
!	Matrix Diagonalizing			   !
!--------------------------------------------------!

  subroutine diag(A,C,S,lambda)
    real(dp) :: A(N,N),C(N,N),S(N,N)
    real(dp) :: lambda(N),work(5*N-1)
    integer :: info
    C = A
    call dsyev('V','U',N,C,N,lambda,work,5*N-1,info)
!    call hegv(C,S,lambda)
    if(info.ne.0) then
      write(*,*)'Diag failed.'
    endif
  end subroutine diag

!--------------------------------------------------!
!	Conmutator of two matrix		   !
!--------------------------------------------------!

 complex(dp) function comm(A,B)
    complex(dp) :: A(N,N), B(N,N)
    dimension :: comm(N,N)
    integer  :: i,j,k
    do i=1,N
       do j=1,N
          comm(i,j) = (0.0_dp,0.0_dp)
          do k=1,N
             comm(i,j) = comm(i,j) + A(i,k)*B(k,j) - B(i,k)*A(k,j)
          enddo
       enddo
    enddo
  end function comm
  
!--------------------------------------------------!
!	Matrix Invertion			   !
!--------------------------------------------------!

   subroutine Inv(A)
    complex(dp) :: A(N,N),C(N,N),D(N,N)
    integer :: i,j
    integer :: info, ipiv(N)
    complex(dp) :: work(5*N-1)       
    C = Id()

  
    call zgesv( N, N, A, N, ipiv,C, N, info )
!   call zgesv( n, nrhs, a, lda, ipiv, b, ldb, info )
    if(info.ne.0) then
      write(*,*)'Diag failed.',info
    endif
    A=C
  end subroutine Inv
  
!--------------------------------------------------!
!	Matrix trace				   !
!--------------------------------------------------!

   complex(dp) function trace(A)
    complex(dp) :: A(N,N)
    integer :: i
    intent(in) :: A 
    trace=0.0_dp
    do i=1,N
      trace=trace + A(i,i)
    enddo    
  end function trace
 
 
end module matrixmod
